#Part 1 The purpose of this notebook is to determine the selective benefit of glucosinolate and flavonoid compounds through plant competition by creating selection gradients. Selection gradients will involve final body mass as the proxy for fitness, but this will be replaced with fitness once the measurement is in. Concentration will be on the x-axis. This analysis will be done in each treatment seperately and account for family and greenhouse location.

#Part 2 The second purpose is to determine if glucosinolates and flavonoids influence suscpetibility to pathogens and if this susceptibility results in increased fitness

Read in and prep data

library(ggplot2)
library(dplyr)
library(lme4)

#read in data
rm(list=ls())
dat<-read.csv("DataSynthesis.csv")

#Remove maple data from the data set for the time being.
dat<-dat[!grepl("aple",dat$Tag, fixed =T),]
#(including maple controls): 
dat<-dat[!dat$treatment=="mcnt",]

#Removing NA's
dat<-dat[!is.na(dat$Sample),]
dat<-dat[!is.na(dat$gluc_Conc),]

#Assign family column
prefamily<-gsub("*.\\|","",dat$Tag)
dat$Family<-gsub("\\-.*","",prefamily)

#remove those with fertilizer treatment, and extra genotypes that are only in the alone treatment. 
dat<-dat[!grepl("i",dat$Sample,fixed=T),]

#Correlation between Glucosinolate and Flavonol Concentration

#Average Duplicates. 
dat2<-dat2 %>% select(-X) %>% group_by(Tag) %>% summarize(ChlorA=mean(ChlorA),ChlorB=mean(ChlorB),gluc_Conc=mean(gluc_Conc),flav_Conc=mean(flav_Conc),Family=first(Family),GA3=first(GA3),treatment=first(treatment),gh_row=first(gh_row),gh_bench=first(gh_bench),GM_TotalLeaf_Area=first(GM_TotalLeaf_Area),comp_number=first(comp_number),ThripsDam=mean(ThripsDam),WhiteFungDam=mean(WhiteFungDam),BlackPathDam=mean(BlackPathDam),Fern=mean(Fern),gh_col=first(gh_col))

package ‘bindrcpp’ was built under R version 3.4.4

#Standardizing greenhouse location The purpose of this is to get one numeric vector which can be used to determine the location of individuals in the greenhouse. It will summarize the distance to the walls of the greenhouse.

plot(dat2$gluc_Conc~dat2$gh_bench*dat2$gh_row)

#Investigating genetic differences
boxplot(gluc_Conc~Family,data=dat2[dat2$treatment=="a",])

boxplot(gluc_Conc~Family,data=dat2[dat2$treatment=="m",])

boxplot(gluc_Conc~Family,data=dat2[dat2$treatment=="gm",])

boxplot(GM_TotalLeaf_Area~Family,data=dat2[dat2$treatment=="a",])

boxplot(GM_TotalLeaf_Area~Family,data=dat2[dat2$treatment=="m",])

boxplot(GM_TotalLeaf_Area~Family,data=dat2[dat2$treatment=="gm",])

#I will nest as bench/row/collumn and determine from there what effects need to be retained or not via anova. 
dat2$GM_TotalLeaf_Area

  [1]  8044  9224  8622 12689 10807  3755 12641  6650 12189  5368  9140
 [12]  6001  3598 10739 15374  9329  5632 10014  5679  6867  4419  2478
 [23]  8360  4515  8532  7050  8540 11180  4570  8811  4884  3500  7054
 [34]  8992  9971  9950    NA  6598 10305  8444  3690 10584 11211  5415
 [45] 14121  3967  3456  7710 10237  5094  4051  7518  7405  6561  6461
 [56]  8124  2024 11096  2912  2988 13042  8765  9731  8220  3558  6068
 [67]  6162 12912  4050 11320  4234 12786  4682  5457  3883  6210 11662
 [78] 12334  2560 10012  9520  8712  3036  5357  7115  6275  5619  4092
 [89]  7113  4083  1260  3065  2365  7518 16025  5094 12801  8684 10321
[100] 11129  8126 13018  5974 11511 13782 11310  2115  4445 10892  4760
[111]  2067 10861 13223  5121  2512  3665  6481  8329  4608 10404 10144
[122]  1736  3480 10920  6780  5988  5647  3010  4075  5648  6410 11052
[133]  2496  8221  9395  4563  8731  4122  3034  9917 10708  6818 10769
[144]  7671  9315 12836 14115 11327  7871 12862  4874 12576  1970 10122
[155] 11602  9984 11558 10532  9761  2304  6574  4122  6356  3725  6434
[166]  7360  4813  8069  9681  9477  7127 10118  4906  3753  6358  7164
[177] 11068  4830 10546 10942 14747  1178  6461  9896 10000 11566 10757
[188] 10358  5601  5934  6281  2520 15967 11153  3862  3092  4796  8456
[199]  9524  9131  9111 12444  8617  5816  3971  5492  5000  4543  7005
[210]  5641 11188  4732  1800  9803  4549 13342  8548  8320 11793  8308
[221]  8092  4920  9979 14448  5970  3477  2646  1592  7780  3448  8475
[232] 12764  9444  2492 10512  4963 12668 11241  7252 11092  9335  6462
[243] 10907  9341  6700  8025  5926  6316 16542  6942  1452 12441 13308
[254]  3612  9806  6969 11863  8976  6572  8411  9304  4889  1748  4812
[265]  7035 13031  4467 13313  2500  3096  6194   972  2122  6370  9504
[276]  4540  6385 10521  7806 15152  8640 10983  9665  8217  4808 10717
[287]  5849  8445 11872 11024  7857  9671  9456  2304  5696  3488 10572
[298]  5992  5434  9953  6930  3854  2550  1691  7227    NA 10215 14991
[309]  4745  3658  7851 14197  9300  4544 10927  4403  8010  8046  8803
[320]  8923  9874  9572  8305  1841 13313  2377  4794 10299  6984  8240
[331]  8309 12474  8522  7805  4149  9960  7687 11250  9343 10524  4928
[342] 11170  5209  7954  6435  7426 13498  7053  8718  6173  5439  4900
[353] 10695 11803  4950 11456 13132  9718  5993 12159 12111  5502  9821
[364]  6174 17238  1923  4768  2088  6471 11503 10236 10909 14955  9739
[375] 12650  2838  5798  5499  7320 16392  7070    NA    NA  9558  9050
[386]  8730  9496 12535  9540 10788 11344 11489  7651  6553  1295  7059
[397]  9187  9993 10512  7858 17606  6060  5167 10399 10159  5003  9306
[408]    NA  7189 10248  6906  6397 11965  8315 12173  8018 11156  7028
[419]  3306 10267  1800 17244  9232  5616  7809 10319 13392 12559 10927
[430]   918  5112  6342  8499  5472

#Effect of glucosinolates & flavonols on fitness (bodymass for proxy)

#removing those without fitness measurements
dat2<-dat2[!is.na(dat2$GM_TotalLeaf_Area),]

#in the maple treatment
ggplot(dat2[dat2$treatment=="m",])+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc))

#in the garlic mustard treatment
ggplot(dat2[dat2$treatment=="gm",])+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc))

#alone
ggplot(dat2[dat2$treatment=="a",])+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc))

#there does not appear to be any on competition  what so ever, however, maybe trends will appear after accounting for gh location. This looks at whether there is an interaction between gluc conc and treatment, which is needed to infer a benefit of glucosinolate concentration. 
library(lattice)

PlotResult<-groupedData(GM_TotalLeaf_Area~gluc_Conc|Family,data=dat2[dat2$treatment=="m",])
PlotResult<-groupedData(GM_TotalLeaf_Area~gluc_Conc|Family,data=dat2[dat2$treatment=="gm",])

plot(PlotResult)

#Event though there does not appear to be a relationship for any given genotype, when we account for their average differences across treatment and greenhouse bench, a trend comes out.


#Visualizing mean result.

dat3<-dat2 %>% group_by(Family,treatment) %>% summarize_if(is.numeric,mean)

#in the maple treatment
ggplot(dat3[dat3$treatment=="m",])+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc,colour=Family))

#in the garlic mustard treatment
ggplot(dat3[dat3$treatment=="gm",])+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc))

#alone
ggplot(dat3[dat3$treatment=="a",])+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc))

#Remove genotypes that had thier competitor die while in the GM treatment. This needs to be done!!!

table(dat2$comp_number)

 1  2  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 
 1  2  1  2  2  1  1  2  1  1  1  2  2  2  2  1  2  1  2  2  1  1  2  2 
26 27 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 51 
 2  2  1  2  1  2  2  2  2  1  1  2  2  2  1  2  1  1  2  1  1  2  2  2 
52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 71 72 73 74 75 76 
 2  2  2  2  1  2  2  1  1  1  1  2  1  1  2  2  2  2  2  2  1  2  2  2 
78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 
 1  1  2  2  2  2  2  2  2  1  1  1  1  2  2 

dat

dat2

#Based on these results, it seems as though i need to include family to avoid psuedo replication, and i need to include bench, but adding greenhouse row does not seem to increase the predictive power that much so i will leave that out. Now i will compare all variables which i am interested in predicting fitness. with these random effects.

#This is the full model, i expect that the influence of gluc conc and flav conc could vary between treatments because of allelopathy, but i have no reason to think that pathogens could influence fitness differently in different treatments. The same goes with ferns. Competition from ferns would affect fitness equally in all treatments, even though there may be more ferns in certain treatments, but this is not what i am testing wtih this analyis. Let us continue with backwards model selection. 

#First lets ensure we are using the correct random effects.
fitfull0<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family), data=dat2)

fitfull<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench), data=dat2)

fitfull2<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench/gh_row), data=dat2)

fitfull3<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench/gh_col), data=dat2)





#Gh_Col is the best predictor for the data by far, so i will use that as the random effects in the model. 

anova(fitfull0,fitfull) #Evidence to use bench at least.

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
object: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
object:     WhiteFungDam + ThripsDam + Fern + (1 | Family)
..1: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
..1:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)    
object 18 7640.8 7712.9 -3802.4   7604.8                             
..1    19 7601.0 7677.2 -3781.5   7563.0 41.724      1  1.051e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(fitfull,fitfull2) #No evidence to use row. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
object: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
object:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench)
..1: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
..1:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_row)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
object 19 7601.0 7677.2 -3781.5   7563.0                         
..1    20 7600.9 7681.0 -3780.5   7560.9 2.1304      1     0.1444

anova(fitfull,fitfull3) #There is strong evidence to use collumn however. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
object: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
object:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench)
..1: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
..1:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col)
       Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)  
object 19 7601.0 7677.2 -3781.5   7563.0                          
..1    20 7597.9 7678.0 -3779.0   7557.9  5.13      1    0.02352 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Modelling fixed effects. 

fitfull3<-lmer(GM_TotalLeaf_Area~treatment*gluc_Conc*flav_Conc+BlackPathDam+WhiteFungDam+ThripsDam+Fern+(1|Family)+(1|gh_bench/gh_col), data=dat2)

#Removing three way interaction
fit.1<-update(fitfull3, ~.-treatment:gluc_Conc:flav_Conc)
anova(fitfull3,fit.1) #Good to remove

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
..1: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
..1:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
..1:     treatment:gluc_Conc + treatment:flav_Conc + gluc_Conc:flav_Conc
object: GM_TotalLeaf_Area ~ treatment * gluc_Conc * flav_Conc + BlackPathDam + 
object:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
..1    18 7600.9 7673.1 -3782.5   7564.9                           
object 20 7597.9 7678.0 -3779.0   7557.9 7.0377      2    0.02963 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.2<-update(fit.1,~.-gluc_Conc:flav_Conc)
anova(fit.2,fit.1) #Good to remove. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.2: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.2:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.2:     treatment:gluc_Conc + treatment:flav_Conc
fit.1: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.1:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.1:     treatment:gluc_Conc + treatment:flav_Conc + gluc_Conc:flav_Conc
      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
fit.2 17 7601.7 7669.8 -3783.8   7567.7                           
fit.1 18 7600.9 7673.1 -3782.5   7564.9 2.7524      1    0.09711 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.3<-update(fit.2,~.-treatment:flav_Conc)
anova(fit.2,fit.3) #Good to remove. 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.3: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.3:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.3:     treatment:gluc_Conc
fit.2: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.2:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.2:     treatment:gluc_Conc + treatment:flav_Conc
      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
fit.3 15 7598.2 7658.3 -3784.1   7568.2                         
fit.2 17 7601.7 7669.8 -3783.8   7567.7 0.5123      2      0.774

fit.4<-update(fit.3,~.-treatment:gluc_Conc)
anova(fit.4,fit.3)  #That significantly affected the predictive power. Did Not Remove 

refitting model(s) with ML (instead of REML)

Data: dat2
Models:
fit.4: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.4:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col)
fit.3: GM_TotalLeaf_Area ~ treatment + gluc_Conc + flav_Conc + BlackPathDam + 
fit.3:     WhiteFungDam + ThripsDam + Fern + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.3:     treatment:gluc_Conc
      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
fit.4 13 7602.4 7654.5 -3788.2   7576.4                           
fit.3 15 7598.2 7658.3 -3784.1   7568.2 8.1828      2    0.01672 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.4<-update(fit.3,~.-flav_Conc)
#It seems as though flav conc has less sample size and so an anova cannot be done.I will refit the model with the same data set but without flav_Conc to test if it is significant. 
fit.4<-lmer(GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam +  
    ThripsDam + Fern +   treatment:gluc_Conc+ (1 | Family) + (1 | gh_bench/gh_col)  ,data=dat2[!is.na(dat2$flav_Conc),])

fit.3<-lmer(GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam +  
    ThripsDam + Fern +   treatment:gluc_Conc+flav_Conc+ (1 | Family) + (1 | gh_bench/gh_col)  ,data=dat2[!is.na(dat2$flav_Conc),])

anova(fit.3,fit.4) #did not remove. 

refitting model(s) with ML (instead of REML)

Data: dat2[!is.na(dat2$flav_Conc), ]
Models:
..1: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
..1:     ThripsDam + Fern + treatment:gluc_Conc + (1 | Family) + (1 | 
..1:     gh_bench/gh_col)
object: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
object:     ThripsDam + Fern + treatment:gluc_Conc + flav_Conc + (1 | 
object:     Family) + (1 | gh_bench/gh_col)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)  
..1    14 7600.5 7656.6 -3786.3   7572.5                           
object 15 7598.2 7658.3 -3784.1   7568.2 4.3304      1    0.03744 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Flav_Conc is definitely a significant predictor on its own. 

fit.4<-update(fit.3,~.-BlackPathDam)
anova(fit.3,fit.4) #do not remove black path dam.

refitting model(s) with ML (instead of REML)

Data: dat2[!is.na(dat2$flav_Conc), ]
Models:
..1: GM_TotalLeaf_Area ~ treatment + gluc_Conc + WhiteFungDam + ThripsDam + 
..1:     Fern + flav_Conc + (1 | Family) + (1 | gh_bench/gh_col) + 
..1:     treatment:gluc_Conc
object: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
object:     ThripsDam + Fern + treatment:gluc_Conc + flav_Conc + (1 | 
object:     Family) + (1 | gh_bench/gh_col)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)   
..1    14 7603.6 7659.7 -3787.8   7575.6                            
object 15 7598.2 7658.3 -3784.1   7568.2 7.4007      1    0.00652 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

fit.4<-update(fit.3,~.-WhiteFungDam)
anova(fit.3,fit.4)  #removing whitefungdam. 

refitting model(s) with ML (instead of REML)

Data: dat2[!is.na(dat2$flav_Conc), ]
Models:
..1: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + ThripsDam + 
..1:     Fern + flav_Conc + (1 | Family) + (1 | gh_bench/gh_col) + 
..1:     treatment:gluc_Conc
object: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + WhiteFungDam + 
object:     ThripsDam + Fern + treatment:gluc_Conc + flav_Conc + (1 | 
object:     Family) + (1 | gh_bench/gh_col)
       Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
..1    14 7596.4 7652.5 -3784.2   7568.4                         
object 15 7598.2 7658.3 -3784.1   7568.2 0.1653      1     0.6844

fit.5<-update(fit.4,~.-ThripsDam)
anova(fit.5,fit.4)  #removing ThripsDam. 

refitting model(s) with ML (instead of REML)

Data: dat2[!is.na(dat2$flav_Conc), ]
Models:
fit.5: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern + 
fit.5:     flav_Conc + (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
fit.4: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + ThripsDam + 
fit.4:     Fern + flav_Conc + (1 | Family) + (1 | gh_bench/gh_col) + 
fit.4:     treatment:gluc_Conc
      Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
fit.5 13 7595.6 7647.7 -3784.8   7569.6                         
fit.4 14 7596.4 7652.5 -3784.2   7568.4 1.2195      1     0.2695

fit.6<-update(fit.5,~.-Fern)
anova(fit.6,fit.5)  #Cannot Remove Fern. 

refitting model(s) with ML (instead of REML)

Data: dat2[!is.na(dat2$flav_Conc), ]
Models:
fit.6: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + flav_Conc + 
fit.6:     (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
fit.5: GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern + 
fit.5:     flav_Conc + (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
      Df    AIC    BIC  logLik deviance Chisq Chi Df Pr(>Chisq)   
fit.6 12 7600.7 7648.8 -3788.4   7576.7                           
fit.5 13 7595.6 7647.7 -3784.8   7569.6 7.134      1   0.007564 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#fit.5 is therefore the best model
#testing to see if the effects remail using the whole data set, because flav_Conc removed 20 observations due to missing values.  
fit.5.Whole<-lmer(GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern +  
(1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc,data=dat2)

summary(fit.5.Whole)

Linear mixed model fit by REML t-tests use Satterthwaite
  approximations to degrees of freedom [lmerMod]
Formula: 
GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern +  
    (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
   Data: dat2

REML criterion at convergence: 7834.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8150 -0.6034  0.0640  0.5573  3.5020 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept)  358959   599.1  
 Family          (Intercept)  102188   319.7  
 gh_bench        (Intercept) 1905371  1380.4  
 Residual                    6951563  2636.6  
Number of obs: 426, groups:  
gh_col:gh_bench, 28; Family, 23; gh_bench, 5

Fixed effects:
                      Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           13830.62    2107.60   110.96   6.562 1.76e-09 ***
treatmentgm           -8384.71    2498.38   125.28  -3.356  0.00105 ** 
treatmentm            -9186.17    2727.11   127.20  -3.368  0.00100 ** 
gluc_Conc             -4026.29    1991.97   126.40  -2.021  0.04537 *  
BlackPathDam            -80.42      27.65   135.22  -2.909  0.00424 ** 
Fern                   -166.07      53.22   135.94  -3.120  0.00221 ** 
treatmentgm:gluc_Conc  4603.32    2514.14   125.19   1.831  0.06948 .  
treatmentm:gluc_Conc   8356.75    2804.11   127.93   2.980  0.00345 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm glc_Cn BlckPD Fern   trtmntg:_C
treatmentgm -0.761                                                
treatmentm  -0.698  0.592                                         
gluc_Conc   -0.944  0.795   0.729                                 
BlackPathDm -0.078  0.036   0.083   0.036                         
Fern        -0.008  0.010  -0.032   0.001 -0.143                  
trtmntgm:_C  0.745 -0.991  -0.581  -0.787 -0.047 -0.008           
trtmntm:g_C  0.669 -0.567  -0.992  -0.707 -0.083  0.023  0.564    

summary(fit.5)

Linear mixed model fit by REML ['lmerMod']
Formula: 
GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern +  
    flav_Conc + (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
   Data: dat2[!is.na(dat2$flav_Conc), ]

REML criterion at convergence: 7443.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8937 -0.6083  0.0186  0.5802  3.4737 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept)  457532   676.4  
 Family          (Intercept)  151352   389.0  
 gh_bench        (Intercept) 1901212  1378.8  
 Residual                    6846335  2616.6  
Number of obs: 406, groups:  
gh_col:gh_bench, 28; Family, 23; gh_bench, 5

Fixed effects:
                      Estimate Std. Error t value
(Intercept)           12500.50    2173.11   5.752
treatmentgm           -7461.78    2531.34  -2.948
treatmentm            -9027.22    2816.04  -3.206
gluc_Conc             -4727.86    2114.45  -2.236
BlackPathDam            -73.25      28.03  -2.613
Fern                   -143.67      53.77  -2.672
flav_Conc              2289.97    1118.24   2.048
treatmentgm:gluc_Conc  3708.32    2548.62   1.455
treatmentm:gluc_Conc   8252.48    2909.58   2.836

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm glc_Cn BlckPD Fern   flv_Cn
treatmentgm -0.767                                            
treatmentm  -0.688  0.590                                     
gluc_Conc   -0.841  0.741   0.668                             
BlackPathDm -0.092  0.041   0.093  -0.002                     
Fern        -0.016  0.010  -0.049  -0.021 -0.140              
flav_Conc   -0.181  0.085   0.069  -0.281  0.113  0.065       
trtmntgm:_C  0.749 -0.991  -0.579  -0.741 -0.050 -0.007 -0.066
trtmntm:g_C  0.659 -0.565  -0.992  -0.646 -0.093  0.042 -0.071
            trtmntg:_C
treatmentgm           
treatmentm            
gluc_Conc             
BlackPathDm           
Fern                  
flav_Conc             
trtmntgm:_C           
trtmntm:g_C  0.562    

#the results did not change. fit.5 is the best model. 

Perhaps i need to include treatment in the nesting because each family may have different amount of individuals in different treatments, which could bias the fitness estimate of that family (the total estimate).

#Diagnostics #it didnt change the result at all and i think it makes more sense so i will keep it.
fit.5.Whole.tr<-lmer(GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern +  
(1 | Family:treatment) + (1 | gh_bench/gh_col) + treatment:gluc_Conc,data=dat2)

plot(fit.5.Whole.tr)

summary(fit.5.Whole.tr)

Linear mixed model fit by REML t-tests use Satterthwaite
  approximations to degrees of freedom [lmerMod]
Formula: 
GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern +  
    (1 | Family:treatment) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
   Data: dat2

REML criterion at convergence: 7835

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.9387 -0.6079  0.0456  0.5633  3.4591 

Random effects:
 Groups           Name        Variance Std.Dev.
 Family:treatment (Intercept)  162669   403.3  
 gh_col:gh_bench  (Intercept)  378481   615.2  
 gh_bench         (Intercept) 1893989  1376.2  
 Residual                     6886655  2624.2  
Number of obs: 426, groups:  
Family:treatment, 69; gh_col:gh_bench, 28; gh_bench, 5

Fixed effects:
                      Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           13801.70    2107.76   228.10   6.548 3.81e-10 ***
treatmentgm           -8413.20    2506.43   448.30  -3.357 0.000856 ***
treatmentm            -8921.50    2745.11   427.00  -3.250 0.001246 ** 
gluc_Conc             -4009.11    1992.29   452.90  -2.012 0.044779 *  
BlackPathDam            -77.64      27.64   459.70  -2.809 0.005175 ** 
Fern                   -166.69      53.28   457.70  -3.129 0.001867 ** 
treatmentgm:gluc_Conc  4630.86    2519.57   451.70   1.838 0.066725 .  
treatmentm:gluc_Conc   8071.45    2820.97   431.40   2.861 0.004425 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm glc_Cn BlckPD Fern   trtmntg:_C
treatmentgm -0.762                                                
treatmentm  -0.699  0.586                                         
gluc_Conc   -0.944  0.794   0.728                                 
BlackPathDm -0.076  0.034   0.079   0.034                         
Fern        -0.009  0.011  -0.033   0.002 -0.140                  
trtmntgm:_C  0.745 -0.990  -0.575  -0.787 -0.044 -0.009           
trtmntm:g_C  0.669 -0.561  -0.991  -0.707 -0.080  0.024  0.558    

confint(fit.5.Whole)

Computing profile confidence intervals ...

                             2.5 %      97.5 %
.sig01                    17.51414  1077.63029
.sig02                     0.00000   722.19382
.sig03                   582.90291  2871.46257
.sigma                  2438.39069  2812.33052
(Intercept)             9768.13584 17915.80109
treatmentgm           -13242.76168 -3507.00853
treatmentm            -14520.80845 -3816.89927
gluc_Conc              -7902.22126  -132.40434
BlackPathDam            -135.12067   -26.45772
Fern                    -269.88524   -61.53812
treatmentgm:gluc_Conc   -307.47814  9488.63612
treatmentm:gluc_Conc    2827.90704 13857.24758

coef(fit.5.Whole)

$`gh_col:gh_bench`
    (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam      Fern
1:1    14457.00   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
1:2    13883.94   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
1:3    13885.42   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
1:4    13461.37   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
2:1    13704.96   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
2:2    13916.79   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
2:3    14103.85   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
2:4    14354.78   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
2:5    13803.48   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
3:1    14208.50   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
3:2    14156.77   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
3:3    14323.58   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
3:4    14037.65   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
4:1    13865.99   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
4:2    14204.20   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
4:3    14204.21   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
4:4    13790.07   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
4:5    13712.26   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
5:1    13441.68   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
5:2    13698.81   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
5:3    13076.56   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
5:4    13747.51   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
5:5    13841.69   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
6:1    13688.24   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
6:2    13061.71   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
6:3    13459.58   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
6:4    13425.92   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
6:5    13740.94   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
    treatmentgm:gluc_Conc treatmentm:gluc_Conc
1:1              4603.321             8356.747
1:2              4603.321             8356.747
1:3              4603.321             8356.747
1:4              4603.321             8356.747
2:1              4603.321             8356.747
2:2              4603.321             8356.747
2:3              4603.321             8356.747
2:4              4603.321             8356.747
2:5              4603.321             8356.747
3:1              4603.321             8356.747
3:2              4603.321             8356.747
3:3              4603.321             8356.747
3:4              4603.321             8356.747
4:1              4603.321             8356.747
4:2              4603.321             8356.747
4:3              4603.321             8356.747
4:4              4603.321             8356.747
4:5              4603.321             8356.747
5:1              4603.321             8356.747
5:2              4603.321             8356.747
5:3              4603.321             8356.747
5:4              4603.321             8356.747
5:5              4603.321             8356.747
6:1              4603.321             8356.747
6:2              4603.321             8356.747
6:3              4603.321             8356.747
6:4              4603.321             8356.747
6:5              4603.321             8356.747

$Family
        (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam
BWPEM1     13711.54   -8384.714  -9186.173 -4026.294    -80.42112
CBMCK1     13955.53   -8384.714  -9186.173 -4026.294    -80.42112
CRSOSO     14016.32   -8384.714  -9186.173 -4026.294    -80.42112
CWRIC2     13850.12   -8384.714  -9186.173 -4026.294    -80.42112
DVGM       13928.01   -8384.714  -9186.173 -4026.294    -80.42112
EFCC2      13738.55   -8384.714  -9186.173 -4026.294    -80.42112
JARI1      13930.65   -8384.714  -9186.173 -4026.294    -80.42112
JBBLB2     13613.39   -8384.714  -9186.173 -4026.294    -80.42112
JBCHY1     13478.37   -8384.714  -9186.173 -4026.294    -80.42112
JWBOY      14049.13   -8384.714  -9186.173 -4026.294    -80.42112
KVEDG1     13573.12   -8384.714  -9186.173 -4026.294    -80.42112
MAVBEL2    13854.25   -8384.714  -9186.173 -4026.294    -80.42112
MHBUR1     13879.65   -8384.714  -9186.173 -4026.294    -80.42112
MHNAT1     13823.37   -8384.714  -9186.173 -4026.294    -80.42112
MSMID1     13901.99   -8384.714  -9186.173 -4026.294    -80.42112
PDVRT1     13941.31   -8384.714  -9186.173 -4026.294    -80.42112
RULEB1     13871.36   -8384.714  -9186.173 -4026.294    -80.42112
SMAKC1     13826.50   -8384.714  -9186.173 -4026.294    -80.42112
SMITH1     14028.71   -8384.714  -9186.173 -4026.294    -80.42112
VRCAN      13767.34   -8384.714  -9186.173 -4026.294    -80.42112
VRPET2     13836.30   -8384.714  -9186.173 -4026.294    -80.42112
VSGARO1    13775.69   -8384.714  -9186.173 -4026.294    -80.42112
WSSWM3     13753.13   -8384.714  -9186.173 -4026.294    -80.42112
             Fern treatmentgm:gluc_Conc treatmentm:gluc_Conc
BWPEM1  -166.0728              4603.321             8356.747
CBMCK1  -166.0728              4603.321             8356.747
CRSOSO  -166.0728              4603.321             8356.747
CWRIC2  -166.0728              4603.321             8356.747
DVGM    -166.0728              4603.321             8356.747
EFCC2   -166.0728              4603.321             8356.747
JARI1   -166.0728              4603.321             8356.747
JBBLB2  -166.0728              4603.321             8356.747
JBCHY1  -166.0728              4603.321             8356.747
JWBOY   -166.0728              4603.321             8356.747
KVEDG1  -166.0728              4603.321             8356.747
MAVBEL2 -166.0728              4603.321             8356.747
MHBUR1  -166.0728              4603.321             8356.747
MHNAT1  -166.0728              4603.321             8356.747
MSMID1  -166.0728              4603.321             8356.747
PDVRT1  -166.0728              4603.321             8356.747
RULEB1  -166.0728              4603.321             8356.747
SMAKC1  -166.0728              4603.321             8356.747
SMITH1  -166.0728              4603.321             8356.747
VRCAN   -166.0728              4603.321             8356.747
VRPET2  -166.0728              4603.321             8356.747
VSGARO1 -166.0728              4603.321             8356.747
WSSWM3  -166.0728              4603.321             8356.747

$gh_bench
  (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam      Fern
1    15861.65   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
2    13504.05   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
3    14199.31   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
4    12947.15   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
5    12640.96   -8384.714  -9186.173 -4026.294    -80.42112 -166.0728
  treatmentgm:gluc_Conc treatmentm:gluc_Conc
1              4603.321             8356.747
2              4603.321             8356.747
3              4603.321             8356.747
4              4603.321             8356.747
5              4603.321             8356.747

attr(,"class")
[1] "coef.mer"

#no heteroscedasticity
plot(fit.5.Whole)

#fairly normal. 
qqnorm(resid(fit.5.Whole))

#Flavonoid is significant, but just barely, there are also alot of samples missing there and I know there is a decent amount of sampling variance there, so i think i will ignore the small positive eeffect of flavonoids on the reduced dataset and continue with the whole one. 

#Conclusion, there is a benefit to producing glucosniolates in the maple treamtment and it is detrimental in the alone treatment. In the garlic mustard treatment, there were mixed results.

#Visualization

library(ggplot2)

coef(fit.5.Whole.tr)

$`Family:treatment`
           (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam
BWPEM1:a      13895.17     -8413.2    -8921.5 -4009.115    -77.64257
BWPEM1:gm     13618.16     -8413.2    -8921.5 -4009.115    -77.64257
BWPEM1:m      13690.67     -8413.2    -8921.5 -4009.115    -77.64257
CBMCK1:a      13788.42     -8413.2    -8921.5 -4009.115    -77.64257
CBMCK1:gm     13955.57     -8413.2    -8921.5 -4009.115    -77.64257
CBMCK1:m      13878.20     -8413.2    -8921.5 -4009.115    -77.64257
CRSOSO:a      13733.00     -8413.2    -8921.5 -4009.115    -77.64257
CRSOSO:gm     13882.56     -8413.2    -8921.5 -4009.115    -77.64257
CRSOSO:m      14130.23     -8413.2    -8921.5 -4009.115    -77.64257
CWRIC2:a      13779.46     -8413.2    -8921.5 -4009.115    -77.64257
CWRIC2:gm     13735.72     -8413.2    -8921.5 -4009.115    -77.64257
CWRIC2:m      13925.59     -8413.2    -8921.5 -4009.115    -77.64257
DVGM:a        13684.40     -8413.2    -8921.5 -4009.115    -77.64257
DVGM:gm       13859.63     -8413.2    -8921.5 -4009.115    -77.64257
DVGM:m        14039.02     -8413.2    -8921.5 -4009.115    -77.64257
EFCC2:a       13707.57     -8413.2    -8921.5 -4009.115    -77.64257
EFCC2:gm      13740.23     -8413.2    -8921.5 -4009.115    -77.64257
EFCC2:m       13794.89     -8413.2    -8921.5 -4009.115    -77.64257
JARI1:a       13907.50     -8413.2    -8921.5 -4009.115    -77.64257
JARI1:gm      13850.78     -8413.2    -8921.5 -4009.115    -77.64257
JARI1:m       13826.83     -8413.2    -8921.5 -4009.115    -77.64257
JBBLB2:a      13827.29     -8413.2    -8921.5 -4009.115    -77.64257
JBBLB2:gm     13763.45     -8413.2    -8921.5 -4009.115    -77.64257
JBBLB2:m      13418.74     -8413.2    -8921.5 -4009.115    -77.64257
JBCHY1:a      13773.18     -8413.2    -8921.5 -4009.115    -77.64257
JBCHY1:gm     13547.51     -8413.2    -8921.5 -4009.115    -77.64257
JBCHY1:m      13444.65     -8413.2    -8921.5 -4009.115    -77.64257
JWBOY:a       13881.99     -8413.2    -8921.5 -4009.115    -77.64257
JWBOY:gm      14127.63     -8413.2    -8921.5 -4009.115    -77.64257
JWBOY:m       13783.77     -8413.2    -8921.5 -4009.115    -77.64257
KVEDG1:a      13591.14     -8413.2    -8921.5 -4009.115    -77.64257
KVEDG1:gm     13670.93     -8413.2    -8921.5 -4009.115    -77.64257
KVEDG1:m      13673.42     -8413.2    -8921.5 -4009.115    -77.64257
MAVBEL2:a     13807.00     -8413.2    -8921.5 -4009.115    -77.64257
MAVBEL2:gm    13801.01     -8413.2    -8921.5 -4009.115    -77.64257
MAVBEL2:m     13846.82     -8413.2    -8921.5 -4009.115    -77.64257
MHBUR1:a      13989.12     -8413.2    -8921.5 -4009.115    -77.64257
MHBUR1:gm     13753.58     -8413.2    -8921.5 -4009.115    -77.64257
MHBUR1:m      13748.85     -8413.2    -8921.5 -4009.115    -77.64257
MHNAT1:a      13812.01     -8413.2    -8921.5 -4009.115    -77.64257
MHNAT1:gm     13744.95     -8413.2    -8921.5 -4009.115    -77.64257
MHNAT1:m      13835.96     -8413.2    -8921.5 -4009.115    -77.64257
MSMID1:a      13896.21     -8413.2    -8921.5 -4009.115    -77.64257
MSMID1:gm     13829.53     -8413.2    -8921.5 -4009.115    -77.64257
MSMID1:m      13805.22     -8413.2    -8921.5 -4009.115    -77.64257
PDVRT1:a      13730.17     -8413.2    -8921.5 -4009.115    -77.64257
PDVRT1:gm     13945.93     -8413.2    -8921.5 -4009.115    -77.64257
PDVRT1:m      13922.60     -8413.2    -8921.5 -4009.115    -77.64257
RULEB1:a      13873.96     -8413.2    -8921.5 -4009.115    -77.64257
RULEB1:gm     13751.97     -8413.2    -8921.5 -4009.115    -77.64257
RULEB1:m      13850.02     -8413.2    -8921.5 -4009.115    -77.64257
SMAKC1:a      13726.88     -8413.2    -8921.5 -4009.115    -77.64257
SMAKC1:gm     13725.68     -8413.2    -8921.5 -4009.115    -77.64257
SMAKC1:m      13946.78     -8413.2    -8921.5 -4009.115    -77.64257
SMITH1:a      13669.61     -8413.2    -8921.5 -4009.115    -77.64257
SMITH1:gm     14089.80     -8413.2    -8921.5 -4009.115    -77.64257
SMITH1:m      13997.47     -8413.2    -8921.5 -4009.115    -77.64257
VRCAN:a       13813.62     -8413.2    -8921.5 -4009.115    -77.64257
VRCAN:gm      13783.00     -8413.2    -8921.5 -4009.115    -77.64257
VRCAN:m       13692.68     -8413.2    -8921.5 -4009.115    -77.64257
VRPET2:a      13874.44     -8413.2    -8921.5 -4009.115    -77.64257
VRPET2:gm     13859.00     -8413.2    -8921.5 -4009.115    -77.64257
VRPET2:m      13685.91     -8413.2    -8921.5 -4009.115    -77.64257
VSGARO1:a     13841.03     -8413.2    -8921.5 -4009.115    -77.64257
VSGARO1:gm    13857.35     -8413.2    -8921.5 -4009.115    -77.64257
VSGARO1:m     13612.91     -8413.2    -8921.5 -4009.115    -77.64257
WSSWM3:a      13836.05     -8413.2    -8921.5 -4009.115    -77.64257
WSSWM3:gm     13545.27     -8413.2    -8921.5 -4009.115    -77.64257
WSSWM3:m      13887.97     -8413.2    -8921.5 -4009.115    -77.64257
                Fern treatmentgm:gluc_Conc treatmentm:gluc_Conc
BWPEM1:a   -166.6875              4630.863             8071.455
BWPEM1:gm  -166.6875              4630.863             8071.455
BWPEM1:m   -166.6875              4630.863             8071.455
CBMCK1:a   -166.6875              4630.863             8071.455
CBMCK1:gm  -166.6875              4630.863             8071.455
CBMCK1:m   -166.6875              4630.863             8071.455
CRSOSO:a   -166.6875              4630.863             8071.455
CRSOSO:gm  -166.6875              4630.863             8071.455
CRSOSO:m   -166.6875              4630.863             8071.455
CWRIC2:a   -166.6875              4630.863             8071.455
CWRIC2:gm  -166.6875              4630.863             8071.455
CWRIC2:m   -166.6875              4630.863             8071.455
DVGM:a     -166.6875              4630.863             8071.455
DVGM:gm    -166.6875              4630.863             8071.455
DVGM:m     -166.6875              4630.863             8071.455
EFCC2:a    -166.6875              4630.863             8071.455
EFCC2:gm   -166.6875              4630.863             8071.455
EFCC2:m    -166.6875              4630.863             8071.455
JARI1:a    -166.6875              4630.863             8071.455
JARI1:gm   -166.6875              4630.863             8071.455
JARI1:m    -166.6875              4630.863             8071.455
JBBLB2:a   -166.6875              4630.863             8071.455
JBBLB2:gm  -166.6875              4630.863             8071.455
JBBLB2:m   -166.6875              4630.863             8071.455
JBCHY1:a   -166.6875              4630.863             8071.455
JBCHY1:gm  -166.6875              4630.863             8071.455
JBCHY1:m   -166.6875              4630.863             8071.455
JWBOY:a    -166.6875              4630.863             8071.455
JWBOY:gm   -166.6875              4630.863             8071.455
JWBOY:m    -166.6875              4630.863             8071.455
KVEDG1:a   -166.6875              4630.863             8071.455
KVEDG1:gm  -166.6875              4630.863             8071.455
KVEDG1:m   -166.6875              4630.863             8071.455
MAVBEL2:a  -166.6875              4630.863             8071.455
MAVBEL2:gm -166.6875              4630.863             8071.455
MAVBEL2:m  -166.6875              4630.863             8071.455
MHBUR1:a   -166.6875              4630.863             8071.455
MHBUR1:gm  -166.6875              4630.863             8071.455
MHBUR1:m   -166.6875              4630.863             8071.455
MHNAT1:a   -166.6875              4630.863             8071.455
MHNAT1:gm  -166.6875              4630.863             8071.455
MHNAT1:m   -166.6875              4630.863             8071.455
MSMID1:a   -166.6875              4630.863             8071.455
MSMID1:gm  -166.6875              4630.863             8071.455
MSMID1:m   -166.6875              4630.863             8071.455
PDVRT1:a   -166.6875              4630.863             8071.455
PDVRT1:gm  -166.6875              4630.863             8071.455
PDVRT1:m   -166.6875              4630.863             8071.455
RULEB1:a   -166.6875              4630.863             8071.455
RULEB1:gm  -166.6875              4630.863             8071.455
RULEB1:m   -166.6875              4630.863             8071.455
SMAKC1:a   -166.6875              4630.863             8071.455
SMAKC1:gm  -166.6875              4630.863             8071.455
SMAKC1:m   -166.6875              4630.863             8071.455
SMITH1:a   -166.6875              4630.863             8071.455
SMITH1:gm  -166.6875              4630.863             8071.455
SMITH1:m   -166.6875              4630.863             8071.455
VRCAN:a    -166.6875              4630.863             8071.455
VRCAN:gm   -166.6875              4630.863             8071.455
VRCAN:m    -166.6875              4630.863             8071.455
VRPET2:a   -166.6875              4630.863             8071.455
VRPET2:gm  -166.6875              4630.863             8071.455
VRPET2:m   -166.6875              4630.863             8071.455
VSGARO1:a  -166.6875              4630.863             8071.455
VSGARO1:gm -166.6875              4630.863             8071.455
VSGARO1:m  -166.6875              4630.863             8071.455
WSSWM3:a   -166.6875              4630.863             8071.455
WSSWM3:gm  -166.6875              4630.863             8071.455
WSSWM3:m   -166.6875              4630.863             8071.455

$`gh_col:gh_bench`
    (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam      Fern
1:1    14467.92     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
1:2    13849.13     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
1:3    13862.44     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
1:4    13439.49     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
2:1    13683.40     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
2:2    13879.28     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
2:3    14072.10     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
2:4    14358.82     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
2:5    13763.75     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
3:1    14201.94     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
3:2    14162.16     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
3:3    14323.11     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
3:4    14002.27     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
4:1    13844.55     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
4:2    14190.87     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
4:3    14181.11     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
4:4    13753.79     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
4:5    13678.42     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
5:1    13379.44     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
5:2    13649.49     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
5:3    13031.63     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
5:4    13701.36     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
5:5    13811.83     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
6:1    13639.40     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
6:2    13012.93     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
6:3    13409.28     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
6:4    13379.28     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
6:5    13718.56     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
    treatmentgm:gluc_Conc treatmentm:gluc_Conc
1:1              4630.863             8071.455
1:2              4630.863             8071.455
1:3              4630.863             8071.455
1:4              4630.863             8071.455
2:1              4630.863             8071.455
2:2              4630.863             8071.455
2:3              4630.863             8071.455
2:4              4630.863             8071.455
2:5              4630.863             8071.455
3:1              4630.863             8071.455
3:2              4630.863             8071.455
3:3              4630.863             8071.455
3:4              4630.863             8071.455
4:1              4630.863             8071.455
4:2              4630.863             8071.455
4:3              4630.863             8071.455
4:4              4630.863             8071.455
4:5              4630.863             8071.455
5:1              4630.863             8071.455
5:2              4630.863             8071.455
5:3              4630.863             8071.455
5:4              4630.863             8071.455
5:5              4630.863             8071.455
6:1              4630.863             8071.455
6:2              4630.863             8071.455
6:3              4630.863             8071.455
6:4              4630.863             8071.455
6:5              4630.863             8071.455

$gh_bench
  (Intercept) treatmentgm treatmentm gluc_Conc BlackPathDam      Fern
1    15835.50     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
2    13469.57     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
3    14149.24     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
4    12924.83     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
5    12629.37     -8413.2    -8921.5 -4009.115    -77.64257 -166.6875
  treatmentgm:gluc_Conc treatmentm:gluc_Conc
1              4630.863             8071.455
2              4630.863             8071.455
3              4630.863             8071.455
4              4630.863             8071.455
5              4630.863             8071.455

attr(,"class")
[1] "coef.mer"

summary(fit.5.Whole)

Linear mixed model fit by REML t-tests use Satterthwaite
  approximations to degrees of freedom [lmerMod]
Formula: 
GM_TotalLeaf_Area ~ treatment + gluc_Conc + BlackPathDam + Fern +  
    (1 | Family) + (1 | gh_bench/gh_col) + treatment:gluc_Conc
   Data: dat2

REML criterion at convergence: 7834.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.8150 -0.6034  0.0640  0.5573  3.5020 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept)  358959   599.1  
 Family          (Intercept)  102188   319.7  
 gh_bench        (Intercept) 1905371  1380.4  
 Residual                    6951563  2636.6  
Number of obs: 426, groups:  
gh_col:gh_bench, 28; Family, 23; gh_bench, 5

Fixed effects:
                      Estimate Std. Error       df t value Pr(>|t|)    
(Intercept)           13830.62    2107.60   110.96   6.562 1.76e-09 ***
treatmentgm           -8384.71    2498.38   125.28  -3.356  0.00105 ** 
treatmentm            -9186.17    2727.11   127.20  -3.368  0.00100 ** 
gluc_Conc             -4026.29    1991.97   126.40  -2.021  0.04537 *  
BlackPathDam            -80.42      27.65   135.22  -2.909  0.00424 ** 
Fern                   -166.07      53.22   135.94  -3.120  0.00221 ** 
treatmentgm:gluc_Conc  4603.32    2514.14   125.19   1.831  0.06948 .  
treatmentm:gluc_Conc   8356.75    2804.11   127.93   2.980  0.00345 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm glc_Cn BlckPD Fern   trtmntg:_C
treatmentgm -0.761                                                
treatmentm  -0.698  0.592                                         
gluc_Conc   -0.944  0.795   0.729                                 
BlackPathDm -0.078  0.036   0.083   0.036                         
Fern        -0.008  0.010  -0.032   0.001 -0.143                  
trtmntgm:_C  0.745 -0.991  -0.581  -0.787 -0.047 -0.008           
trtmntm:g_C  0.669 -0.567  -0.992  -0.707 -0.083  0.023  0.564    

ggplot(dat3)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=gluc_Conc,colour=treatment))+
  geom_abline(intercept=13830.62,slope=-4026.29,colour="red")+
    geom_abline(intercept=13830.62-9186.17,slope=-4026.29+8356.75,colour="blue")

ggplot(dat3)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=BlackPathDam,colour=treatment))

ggplot(dat2)+
  geom_point(aes(y=GM_TotalLeaf_Area,x=Fern))+
  geom_abline(intercept=13830.62,slope=-166.07)

We wouldnt really expect to have the power to detect a GXE effect on fitness here, because each genotype will have a specific degree of gluc concentration. I dont think that i can do a model using the variables from leaf data, because the dependent variable is the same for each leaf, they come from the same plant.

#——————————————————————- GROWTHRATE #——————————————————————-

#Repeat this analysis on RGR1

GrowthDatMean<-left_join(dat3,GrowthDatMean,by=c("Family","treatment"))

Column `Family` joining character vector and factor, coercing into character vectorColumn `treatment` joining factors with different levels, coercing to character vector

#Modelling relative growth rate 1. Pathogens and fern abundance was removed from this analysis because they did not appear this early in the experiment.

#Modelling fixed effects. Ignoring flavonoids for the time being. 

fitfull3<-lmer(RGR1~treatment*gluc_Conc+(1|Family:treatment)+(1|gh_bench/gh_col), data=GrowthDat)


#Removing three way interaction
fit.1<-update(fitfull3, ~.-gluc_Conc:treatment)
anova(fitfull3,fit.1) #Good to remove

fit.2<-update(fit.1,~.-gluc_Conc)
anova(fit.2,fit.1) #Good to remove. 

fit.3<-update(fit.2,~.-treatment)
anova(fit.2,fit.3) #Significant

Treament was the only predictor for growth rate 1.

#Modelling on growth rate 2.

anova(fit.2,fit.3) #Significant

refitting model(s) with ML (instead of REML)

Data: GrowthDat
Models:
fit.3: RGR2 ~ (1 | Family:treatment) + (1 | gh_bench/gh_col)
fit.2: RGR2 ~ treatment + (1 | Family:treatment) + (1 | gh_bench/gh_col)
      Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)    
fit.3  5 -11315 -11295 5662.6   -11325                             
fit.2  7 -11331 -11303 5672.7   -11345 20.214      2  4.079e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Treament was the only predictor for growth rate 2.

#Modelling on growth rate 3.

anova(fit.6,fit.7) #Not Significant

refitting model(s) with ML (instead of REML)

Data: GrowthDat[!is.na(GrowthDat$ThripsDam), ]
Models:
fit.7: RGR3 ~ (1 | Family:treatment) + (1 | gh_bench/gh_col)
fit.6: RGR3 ~ Fern + (1 | Family:treatment) + (1 | gh_bench/gh_col)
      Df    AIC    BIC logLik deviance  Chisq Chi Df Pr(>Chisq)  
fit.7  5 -11574 -11554 5792.0   -11584                           
fit.6  6 -11575 -11551 5793.6   -11587 3.1064      1    0.07799 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Modelling on growth rate 4. This one include the other variables as they were now involved in the experiement around this time.

summary(fit.2)

Linear mixed model fit by REML ['lmerMod']
Formula: 
RGR4 ~ treatment * gluc_Conc + Fern + BlackPathDam + (1 | Family:treatment)
   Data: GrowthDat

REML criterion at convergence: -11683

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-4.0620 -0.4092  0.0438  0.5709  5.4749 

Random effects:
 Groups           Name        Variance  Std.Dev. 
 Family:treatment (Intercept) 2.418e-16 1.555e-08
 Residual                     1.090e-14 1.044e-07
Number of obs: 408, groups:  Family:treatment, 69

Fixed effects:
                        Estimate Std. Error t value
(Intercept)           -1.595e-07  8.234e-08  -1.937
treatmentgm            2.462e-07  1.021e-07   2.410
treatmentm             2.626e-07  1.111e-07   2.363
gluc_Conc              1.715e-07  8.105e-08   2.116
Fern                   4.398e-09  1.991e-09   2.209
BlackPathDam          -5.202e-09  1.105e-09  -4.709
treatmentgm:gluc_Conc -2.345e-07  1.026e-07  -2.285
treatmentm:gluc_Conc  -3.034e-07  1.138e-07  -2.665

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm glc_Cn Fern   BlckPD trtmntg:_C
treatmentgm -0.804                                                
treatmentm  -0.744  0.597                                         
gluc_Conc   -0.992  0.799   0.736                                 
Fern        -0.005  0.000  -0.046  -0.005                         
BlackPathDm -0.078  0.035   0.100   0.034 -0.102                  
trtmntgm:_C  0.785 -0.991  -0.584  -0.790  0.006 -0.048           
trtmntm:g_C  0.712 -0.571  -0.992  -0.715  0.037 -0.100  0.567    

#The next step to support this analysis would be to determine if glucosinolate content explains competitive ability against maples.

#Look at if there is a negative correlation between biomass in the maple and alone treatment.

#—————————————– #Part 2

#assessing

#Significant testing (doing seperately because flav conc and gluc conc are correlated)

#WhitePathDam

#BlackPathDam

summary(fit.4)

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: BlackPathDam ~ flav_Conc + (1 | Family/Tag) + (1 | gh_bench)
   Data: dat[!is.na(dat$flav_Conc), ]

     AIC      BIC   logLik deviance df.resid 
  3504.2   3526.8  -1747.1   3494.2      672 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0421 -0.4065 -0.2289  0.0876  4.0742 

Random effects:
 Groups     Name        Variance Std.Dev.
 Tag:Family (Intercept) 0.971171 0.98548 
 Family     (Intercept) 0.052429 0.22897 
 gh_bench   (Intercept) 0.003662 0.06051 
Number of obs: 677, groups:  Tag:Family, 423; Family, 23; gh_bench, 5

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.6816     0.1579  10.651  < 2e-16 ***
flav_Conc    -1.1150     0.1644  -6.781 1.19e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr)
flav_Conc -0.856

#ThripsDam

summary(fit.4) #Flav conc is again, a very significant predictor. 

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: BlackPathDam ~ flav_Conc + (1 | Family/Tag) + (1 | gh_bench)
   Data: dat[!is.na(dat$flav_Conc), ]

     AIC      BIC   logLik deviance df.resid 
  3504.2   3526.8  -1747.1   3494.2      672 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0421 -0.4065 -0.2289  0.0876  4.0742 

Random effects:
 Groups     Name        Variance Std.Dev.
 Tag:Family (Intercept) 0.971171 0.98548 
 Family     (Intercept) 0.052429 0.22897 
 gh_bench   (Intercept) 0.003662 0.06051 
Number of obs: 677, groups:  Tag:Family, 423; Family, 23; gh_bench, 5

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)   1.6816     0.1579  10.651  < 2e-16 ***
flav_Conc    -1.1150     0.1644  -6.781 1.19e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
          (Intr)
flav_Conc -0.856

summary(fit.1) #Flav conc and treatment predict fern abudnance. 

Generalized linear mixed model fit by maximum likelihood (Laplace
  Approximation) [glmerMod]
 Family: poisson  ( log )
Formula: Fern ~ treatment + flav_Conc + (1 | Family) + (1 | gh_bench/gh_col)
   Data: dat[!is.na(dat$flav_Conc), ]

     AIC      BIC   logLik deviance df.resid 
  1804.4   1836.1   -895.2   1790.4      680 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.5073 -0.6085 -0.2213 -0.1154 16.0269 

Random effects:
 Groups          Name        Variance Std.Dev.
 gh_col:gh_bench (Intercept) 3.636    1.907   
 Family          (Intercept) 1.701    1.304   
 gh_bench        (Intercept) 1.153    1.074   
Number of obs: 687, groups:  gh_col:gh_bench, 28; Family, 23; gh_bench, 5

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -2.1683     0.7522  -2.883 0.003944 ** 
treatmentgm  -0.2786     0.1498  -1.859 0.063030 .  
treatmentm    0.4068     0.1122   3.625 0.000289 ***
flav_Conc    -0.6571     0.2737  -2.401 0.016365 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
            (Intr) trtmntg trtmntm
treatmentgm -0.148                
treatmentm  -0.165  0.500         
flav_Conc   -0.325  0.221   0.252